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Abstract 

Random linear systems over the Galois Field modulo 2 have an 
interest in connection with problems ranging from computational op- 
timization to complex networks. They are often approached using 
random matrices with Poisson-distributed or finite column/row-sums. 
This technical note considers the typical rank of random matrices 
belonging to a specific ensemble wich has genuinely power-law dis- 
tributed column-sums. For this ensemble, we find a formula for calcu- 
lating the typical rank in the limit of large matrices as a function of the 
power-law exponent and the shape of the matrix, and characterize its 
behavior through "phase diagrams" with varying model parameters. 



1 Introduction 

This technical note presents the calculation of the typical rank of Boolean 
random matrices with power-law distributed column-sums. The specificity 
of this calculation is that it applies to genuinely power-law matrices, without 
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finite cutoffs in the distribution. Before presenting the results, we will give 
a brief description of the context that motivates the calculation. 

Random matrices with Boolean entries are often simple to treat, which 
makes them important in many paradigmatic problems of different branches 
of science. For example, in computer science, they define the so-called ran- 
dom XOR-SAT problem [TH3] , the simplest of an important class of opti- 
mization problems at the interface of statistical physics [HE] and computer 
science [6HH]. The XOR-SAT problem consists in finding a solution to the 
set of linear equations of N Boolean variables and M equations Aa = f over 
the Galois Field of order 2 (usually indicated as OF (2)), where the matrix A 
is extracted from a prescribed ensemble of Boolean matrices. 

The typical properties of the linear systems can be computed in the limit 
of large matrices and fixed density of constraints 7 = M/N. For random ma- 
trices with constant row-sums (and thus Poisson-distributed column-sums), 
the "order parameter" 7 plays a crucial role for the solution space of the cor- 
responding random XOR-SAT problem [I]. With increasing 7, the random 
XOR-SAT presents three different regimes with some features of a thermody- 
namics phase [9]. For 7 < 7^ a solution can be typically found by removing 
iteratively all variables present in only one equation (trivial pivots in the 
language of Gaussian elimination [UEEU]). m this case, it can be shown that 
the solution space is composed of only one cluster. For 7^ < 7 < j c matrices 
have typically a non-empty "core" (the remaining part of the matrix after 
the recursive elimination of the trivial pivots) and finding a solution requires 
a number of iterations proportional to the cube of the size of the core [10 ]. 
Here, the solution space is split into many well separated clusters. Finally, 
for 7 > 7 C in the typical case solutions cannot be found (i.e. the solution 
space is empty). 

In the field of complex networks, Boolean matrices are used to represent 
empirical systems with many interacting agents: each agent is labelled with 
an integer and the entry of the matrix Aij is equal to one only if agent i 
interacts with agent j, and zero otherwise. For instance, properties of the 
matrix A are useful to control graph properties like hyperloops or critical 
sets of independent nodes [IT] . In order to study the typical properties of 
such a system, it is necessary to define an ensemble of matrices which con- 
serves characteristic properties of the empirical case. Of particular interest 
are matrices with a power-law distribution of column-sums, which are typical 
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of many empirical graphs [T21fi3] . 

We have previously introduced a simple and analytically treatable Boolean 
random matrix ensemble with a power-law distribution of the column-sums 
p(k) ~ k~P and tunable HSHH]. This paper describes an analytical ap- 
proach to the problem of the typical rank over GF(2) of random matrices 
belonging to the this ensemble and compares the results to a numerical eval- 
uation. Previous approaches of this kind were applied to similar and more 
sophisticated models, but were limited to distributions of the row/column- 
sums with Poisson |3] or regular tails [18], or with power-law tails with a 
finite cut-off PUCES]. 

The calculation presented here is similar to the replica calculation for 
spin-glasses [20J. It allows to find a formula for the typical rank in the limit 
of large matrices as a function of the model parameters 7 and 0, which allows 
to derive interesting phase diagrams. In particular, we estimate a second 
order transition in the typical rank varying the parameter 7. We compares 
the results with the structure of solution space obtained numerically. These 
results are resumed by interesting phase diagram for the behavior of the 
linear system with varying density of constraint 7 and power-law exponent 
0. 

2 Matrix Ensemble 

This paragraph briefly describes the matrix ensemble. A more exhaustive 
characterization can be found in p^6|[T7] . 

The matrix ensemble (Fig. [TJ was originally formulated as a null model 
for (biological) transcriptional regulatory networks. It is defined by the fol- 
lowing generative algorithm. For each column of A, (i) throw a bias from a 
prescribed probability distribution 7tM{d6) and (ii) set the column elements 
of A to be or 1 according to the toss of a coin with bias 9. Since each 
column is thrown independently, the resulting probability law is 

N 1 

p(A) = n / e? r=i Ai] (i - ^) Ef=i(i_A ' ) (i) 

Note that only columns are independent, while the row elements are not 
independent, but symmetric by permutations. 
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Figure 1: Schematic representation of the matrix ensemble. The probability 
that a variable is involved in k constraints is asymptotically proportional 
to a power-law p(k) ~ in the limit of large matrices. Vice versa, the 
probability that a constraint contains s variables is a Poisson distribution 
p(s) ~ ^r e_A 5 where A > is defined in the text. 



To complete the model, one has to specify the choice for ir M (d9), which de- 
termines the behavior of the graph ensemble. To obtain a power-law column- 
sums distribution we choose the two-parameter distribution 

7r M (de) = zrfe-Pxfo ifo, (2) 

where a > and ft > 1 are free parameters, x^°l i] is the characteristic 
function of the interval f-^, l], taking the value one inside the interval and 
zero everywhere else, and Zm = ( M /^^ ~ 1 j s the normalization constant. 
The function of Eq. [2] gives a power-law tail to the column-sums distri- 
bution. Conversely, the cutoff on 9 defined by a poses a constraint on the 
number of nodes with low degree, and will be used to control the probability 
to extract a node with small k. In the limit of large graphs (i.e. in the limit 
M, N — > oo, with M/N = 7 < 00) the probability to extract a matrix with 
hi ones in the i — th column is asymptotically 

°o f0 0. ki -t 

i=i ^° i ' 

where 

VToo(^) = {ft~ l)« /3 ~ 1 X[a,oo)t" /3 dt (4) 

is the limit of the distribution in Eq. [2j Eqs. |3] and H] imply that the probabil- 
ity to have a column with k ones and the probability to have a row with s ones 
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in the limit of the large graphs are respectively p c {k) = ye t TT OD (dt) w 
k~P, and p r (s) = 7fe" A , where A = 7 J °° tn^dt). 

Fig. [2] reports the distribution of the nonzero entries of matrices extracted 
from the ensemble described by Eq. |3l for different values of a and /3. As 
expected, the column-sums (top) follow a power-law distribution while the 
distribution of row-sums (bottom) follow a Poisson distribution. For 1 < 
< 2 the mean row-sum depends on the dimension of the system as /i = 

O^j N 2 ~P ', while for > 2, the mean value of the distribution is 

independent of the size of the system and it is /i — §5? 



3 Calculation of the Typical Rank 

We will now consider the rank of a matrix belonging to the ensemble de- 
scribed in the previous paragraph. There are different methods for comput- 
ing the rank of a given matrix A. Here, we exploit the calculation of the 
number of solutions of the corresponding homogeneous linear system 

where a G {0, 1}^ and 5 (a) , mod 2 % is different from zero only if a = (mod 2). 
Since linear algebra applies, the number of solutions of the homogeneous 
system over the finite field GF(2) can be expressed in terms of the dimension 
of the kernel of matrix A 

AT (A) = 2 null( - 4) . 
Using the rank-nullity theorem 

rank(^l) + null (.4) = N, 

the typical rank of random matrices will be 

(rank(^))=JV-{lDg 2 jV(^)), 

where the average (•) is carried over the matrix ensemble in Eq. [TJ 
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Figure 2: Distribution of nonzero entries for the matrix ensemble (Eq. [3]), at 
(3 = 1.8 (left) and (3 = 2.8 (right). As reported in the text, the column-sums 
(top) follow a distribution having a power-law tail with exponent f3. The 
dashed (green) line is a guide to the eye. On the other hand, the distri- 
bution the row-sums (bottom) follows a Poisson-like distribution with mean 
depending on the value of the parameter /3. 
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In order to calculate the logarithm of the number of solutions we use the 
known limit 



log* = lim^— i (5) 

fe->o k 



where X is a generic random variable [20|[2Tj . In principle, using the above 
limit, it is possible to calculate the average (logAf) knowing the function 
4>{k) = (X k ^, where k is a real parameter. However, the calculation of the 
function <f)(k) for any real k is typically hard. 

As proposed in [20l[2T], a feasible protocol to compute 4>(k) consists in 
calculating the k — th moment of the random variable (X k } (i.e. evaluate 
<f)(k) for integer values of k) and then finding by interpolation a reasonable 
extension for any real k. In many cases it is possible to find a well- 

behaved extension of the function <f)(k), but this is not generally true 



Thus, we are interested in the calculation of the k — th moment of the 
number of solutions 4>{k) = (Af(A) k y As reported in Appendix IA1 we find 



{m s } V 



M 
m 



£6/ £m s ]<snr[ 

Te[k] \Se[k] 



-, N 

, (6) 



where the sum is carried over 2 k integer variables labelled by an element 
of [k] (i.e. the set of all possible subsets of {1,2, ... ,k}) constrained by 
J2se[k] m s — 1 ("replica" indices), and }fl[ is equal to one if the cardi- 
nality of the set Q is odd and zero otherwise. The function £m(^) : = 
J TiM^dO) (1 — 29) h is related to the moments of the column-sum distribu- 
tion of the random matrices extracted from the ensemble in Eq. [2J Note that 
Eq. [6] is not an approximation but it is valid for any M, N < oo. 

In the limit M — > oo at fixed x = h/M the function ^M(h) can be written 

as 

£(x) = lim £(/i/M) = f TToo^e-^, 

M— too J 

where ir^dt) = liniM^oo ^M^dO) (see Eq. HI). It is immediate to observe that 
£ (x) is the moment-generating function of Eq. [1] 



Pc (k) = (-fk\ 



dx k 
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Using the denning expression for the ensemble (Eq. [JJ), Eq. [6] can be 
rewritten as 



y—km 



/ [dx] exp N< 7 S(x < s) 
^ I se[k] 



+ log 



re[fe] \5e[Jfc] 



+ o(l/iV) , (7) 



where the integration is carried over the rescaled variables x$ = mg/N (with 

the constraint J2se[k] x s = 1 ) and Sse[fc] ©fas) = Ese[fc] -zs log 2:5 is the 
Shannon entropy. The above expression diverges exponentially with the di- 
mension N of the matrices, and thus it is possible to use the saddle point 
approximation. In order to compute the saddle point, it is necessary to find 
the maximum of Eq. [7J varying x$, i.e. it is necessary to solve a system of 
a 2 k variables for any integer k. Obviously, this is unfeasible and one must 
impose a symmetry ansatz for the saddle point solution in order to reduce 
the number of variables. 

The simplest hypothesis it that the most symmetric solution would dom- 
inate (in the theory of glassy systems this solution is usually called replica 
symmetric (RS) solution) 



x$ = x 

x s = 1 - (2 k - l)x, 5^0, 

where all variables are equal, except one in order to satisfy the constraint 
J2se[k] x s = 1- Here, the variable x plays the same role of the "Edward- 
Anderson" order parameter in the Spin Glass theory [20]: for x — 0, the 
total entropy ^ Sg [ fc ] &(xs) is exactly zero and then only one state, i.e. the 
most symmetric state, dominates the saddle point in Eq. [71 On the con- 
trary, for x = 1, the total entropy assumes the highest possible value and 
then many different states contribute to the saddle point in Eq. [7J 



Using the RS ansatz, the asymptotic behaviour of the k — th moment of 
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the number of solutions can be written as 



log Aft (A) = _^ bg 2 + ^ r k _ + 
N x I 

+ 7 6(1 -x) + log [1 + (2 fe - l)£(2 k - l x)] }. (8) 

It is important to observe that the variable k in Eq. [S] can assume any real 
value and it can be considered as a possible extension of the Eq. [HI in the limit 
of large matrices. Eq. [8] depends directly on the chosen symmetry ansatz and 
is not guaranteed to be consistent. In our case, we will show that Eq. [8] gives 
results that agree with numerical results. 

We can now take the limit k — >■ 0. Thus we have 

(rank (A)) (log 2 A/o) f ~ , \ t fx 



1 - lim = max | 7 6 (x) - 7 + , 0) 

where ©o(^) = — £ log x + x. The above equation can be used directly to find 
the typical rank of the matrices extracted from the matrix ensemble proposed 
in Eq. [TJ Fig. |3] compares the theoretical prediction of the typical rank with 
simulations. It is possible to observe that, independently of the choice of the 
parameters a and /3, the theoretical prediction is in good agreement with the 
simulations. 

Interestingly, the theoretical prediction of the rank (Eq. |9]) can have a 
second order discontinuity varying the density of constraints 7, due to the fact 
that the value of the RS order parameter x which maximize the expression 
in Eq. |9]can have a jump (Fig. H]). In particular, we find that for any j3 > 2 
there exists a critical value a c (/3) such as for a < a c ((3) there are no jumps 
varying the parameter 7. Instead, for a > a c ((3), it is possible to identify a 
critical value 7 c (/3) in which x has a jump. On the contrary, for 1 < (3 < 2 a 
discontinuity is always present. 

The presence of a second order discontinuity of the typical rank is a signal 
of the fact that the totally symmetric solution (RS solution) is no longer 
valid (even if it may still be a good approximation for the calculation of the 
typical rank) caused by a spontaneous symmetry breaking of the solution 
space in many well-separated clusters [20J. In this less symmetric 

solution (called replica symmetry breaking (RSB) solution) dominates the 
saddle point in Eq. [71 We did not explore analytically this regime. 
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Figure 3: Distribution of the typical rank obtained from simulation with 
N = 500 and /3 = 1.8 (left) or (3 = 2.6 (right), varying the parameter 7. As 
shown in the figures, the numerical data are in agreement with the theoretical 
prediction obtained by Eq. |9j The deviation for small values of 7 is due to 
the small system size. 
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Figure 4: Value of the RS order parameter x that maximizes Eq. [9] at fixed 
p. For < 2 (left), for any value of a there exists a critical value of 7 in 
which the value of x at the maximum has a jump. For > 2 (right) and a 
sufficiently small, the value x max does not have any discontinuity. Otherwise, 
it is possible to identify a 7 C (that depends on a and (3) for which the value 
of x max has a jump. 
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(a) X2 + X'i = 0, mod 2 

(fj) Xi + x'2 + X3 = 0, mod 2 
(7) x i + x 3 + x 4 = 0, mod 2 



Figure 5: Factor graph representation of the XOR-SAT problem. In the 
sketch the variables (columns) are represented by circles and the constraints 
(rows) by rectangles. 

4 Leaf Removal and Organization of the So- 
lution Space 

As described in the previous paragraph, the typical rank of A is related to 
the total number of solutions of the linear system Aa = 0. In particular, we 
found an analytical expression for the typical rank which has sharp transi- 
tions when the parameters that define the matrix ensemble vary continuously. 
As previously discussed, these transitions are related to the clusterization of 
the solution space. This paragraph focuses on the geometrical organization 
of the solution space of the linear system Aa = f (the XOR-SAT prob- 
lem) and the comparison between numerical evaluations and our theoretical 
predictions. A general introduction to this problem can be found in [HllU[ l25j. 

A system of linear equations in GF(2) can be conveniently represented by 
factor graphs, defined by the matrix A, in which variables and constraints 
correspond to distinct types of nodes. If the variable i is present in the 
constraint a, a link (i, a) is drawn in the factor graph (Fig. 

Following [UE3HIE5], ^ is possible to obtain a precise definition of clusters 
of solutions using the so-called "leaf removal" algorithm. The leaf removal 
algorithm is an iterative algorithm used to gradually eliminate all trivially 
constrained variables (called trivial pivots in the language of Gaussian elim- 
ination). It is easy to prove that when a variable (called "leaf") is connected 
to only one constraint, it is always possible to choose its value such that the 
constraint is always satisfied (e.g. variable 4 in Fig. The leaf removal 
algorithm is based on this evidence and it is defined as follows: (i) pick a 
variable that appears only in one constraint (leaf) and (ii) remove it together 
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with the only constraint it is connected to. The process is iterated until no 
leaves remain. The part of the factor graph that cannot be removed by leaf 
removal iteration is called "core" and does not depend on the order in which 
the leafs are removed. In this case, the order parameter of the reduced linear 
system will be 



leave = -T7 , (10) 

-i V r nrp 



i.e. the density of constraints that are not trivially satisfied. 

The presence of the core is related to the clusterization of the solution 
space. If icore = (no core is present), the problem to find a solution of 
the linear system Aa = f is trivial (the complete solution can be found by 
running the leaf removal in reverse direction, in a scheme usually called leaf 
reconstruction) and the solution space is composed of only one cluster. If 
< Jcore < 1, the core is not trivial (but not over-constrained) and each 
solution of the linear system reduced to the core variable defines a single 
cluster. All the solutions built from a core solution by leaf reconstruction be- 
long to the same cluster. Finally, for 7 core > 1 the reduced linear system for 
the core variables is over-constrained, so that no solutions are typically found. 

Fig. E] reports the curves of the typical 7 core varying the density of con- 
straints 7 obtained by numerical simulations of the leaf removal algorithm. 
As predicted in the previous paragraph, the presence of a non over-constrained 
core depends on the choice of the parameter /3. For < 2 (left panel), vary- 
ing the parameter 7 it is always possible to identify three regimes: an empty 
core phase (■y C ore = 0), a non over-constrained core phase (7 core < 1) and an 
over-constrained core phase (jcore > 1)- On the other hand, for (3 > 2 (right 
panel) the not over-constrained (7 core < 1) core is present only for a suffi- 
ciently large. All these results are resumed in the phase diagrams obtained 
from in Fig. [7J 



5 Conclusion 

In conclusion, we have presented a simple calculation of the typical rank 
of random matrices with power-law distributed column-sums on the Galois 
Field of order 2. The matrices can describe a graph or a sparse linear system 
for Boolean variables. The calculation is based on a fairly standard replica- 
like approach, where we compute the generic k-th moment of the number of 
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Figure 6: Numerical simulation of the 7 cor e varying the parameter 7, for 
different value of a and 0. For (3 < 2 (left), it is always possible to find 
the critical value 7 C of inversion of the core. For /3 > 2 (right), only for a 
sufficiently large it is possible to find the critical value j c . In parenthesis the 
theoretical predictions. 



solutions of the associated linear system and we consider the limit k — > of 
its analytical extension in the maximally symmetric case. 

Differently from other models present in the literature [HQIllIIElHH], t ne 
simplicity of the matrix ensemble [16J that we employ here allows to find an 
analytical expression for the typical rank without having to impose any cutoff 
on the power-law distribution. As shown in Figs. [3j the typical rank calcu- 
lated with our method is in fairly good agreement with the numerical results. 
We find that, as usually happens in this kind of models [UEEO] the typical rank 
can have a second order discontinuity with increasing density of constraints 
7. This discontinuity is related to the clusterization of the solution space 
in many well separated clusters of the related XOR-SAT problem [23]. Our 
result indicates that the same phenomenology can exist in presence of truly 
power-law tails. 

More in detail, since the matrix ensemble is defined as a function of the 
model parameters a, which sets a lower cutoff on the row-sums and (3, the 
exponent of the column-sum distribution, one can study the variation of 
this threshold with "phase diagrams" where these parameters vary together 
with the density of constraints. Specifically, the presence of the typical rank 
discontinuity at 7 = 7 C depends on the choice of a and (3. For (3 < 2 the 
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Figure 7: Phase diagrams obtained from simulation with M = 500 and 
iV = M/7. In the plots the results of the simulation are compared with the 
theoretical predictions (solid lines). The top panels contain fixed (3 phase 
diagrams for /3 > 2 and /3 < 2, while the bottom panel is a fixed a phase 
diagram. The largest errors in the phase diagrams arise around the critical 
value a c (top panels) and (3 C (bottom panel). This can be explained observing 
that near such critical values the matrix core is very small compared to the 
finite size fluctuations. 
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discontinuity exists for any choice of a. Otherwise, for /3 > 2, it is possible 
to identify a critical value a c (/3) such that only for a > a c ((3) a critical value 
7 C exists. 

The role played by a in this model is similar to that played by the con- 
straint connectivity K in the i^-XOR-SAT problem [2]. In this case, the 
row-sum of the matrix is equal to K, and the clustering of solution is possi- 
ble only if K > 2. In our case, it is simple to verify that only for < 2 the 
fraction of rows with two nonzero entries always vanishes for every a in the 
large N limit. Thus, we speculate that the density of rows with two or less 
nonzero entries may become important and affect in some cases the phase 
diagram for (3 > 2, causing the observed lack of the clusterization regime. 

Finally, the approach presented here, suitably generalized, may be useful 
to study self-organizing properties of systems with many interacting agents, 
where similar threshold phenomena can emerge as a function of the prop- 
erties of the network that defines the agent interactions. In this case, the 
parameters of the matrix ensemble represent tunable quantitative topolog- 
ical properties of the interaction network such as the connectivity and the 
density of interactions. 

A Calculation of (Afg(A)) 

In this appendix we explicitly calculate the k — th moment of the number 
of solutions of the homogeneous linear system Nq(A) = ^2^6 (Aa), with 
a G {0,1} and 5(a) = 1 only if a — 0. Let p(A) a generic probability 
distribution for the random matrix A: thus the k — th momentum can be 
written as 

M k / N 

XG{0,l} iV ®{0,l}' : ^e{0,l} M ®{0,l} iV j=la=l \i=l 

For simplicity, in the rest of the appendix we use the convention 

i G N, i — 1, . . . , N (position of the row) 
j G N, i = 1, . . . , M (position of the column) 
aGN, 2 = 1, . . . , k (number of the "replica") 
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Use probability distribution for our model (Eq. [T]), the expression of the 
k — th moment will be 



WW = n 



1 -|- ( — l^YUAjiXia 



j, a 



2-"' EE 



A 



JJ(l + (_l)Ei>»iAa) 



where we used the explicit representation of the Kronecker delta for binary 
variables 

At this level, it is possible to exchange the sums over A and the integration 
to obtain 



n 



j2 n i 1 + (-i) EiaiXia ) ( x - ^) i_oi 



The last term does not depend explicitly on j and then we above expression 
can be rewritten as 



-l M 



a£{0,l} N a « 



Now, using the identity 



n c 1 +/(«))= e n/H 

SC[fc] «es 
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where [k] is the set of all the possible subsets of {1, k}, the above 
expression becomes 



En 

SC[k] i 



£ ((-i^^ra-ft) 1- * 

ae{o,i} 



M 



It easy to observe that the last term can be directly calculated. Thus, after 
a sum over a we obtain 

E ((-i) E -^r (i - = i - 2M (V £ xj) , 

S{0,1} V aes J 



cre{0,l} 

where 5(1,0") equals 1 if and only if a — 1, and 



En 

5C[fc] i 



1-2M i,£x, 



ae5 



M 



In order to complete the calculation, it is necessary to expand the term inside 
the curly brackets. Let {m^} the set of 2 fc variables such that J2se[k] m s — 
M. Thus we have 



W4 = 2-'"EE(") 



X {m s } 



n 



ttm(^) n 

5C[fe] 



l-20 i <j(l,£x ia 



where (^) is the multinomial. Using the simple identity 



1-2M (£x ia 



'1 — 2^) 5 (^ aesXia ) m ' s , 



we obtain 
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where we used the notation 

Xj(iS) := ^ Xj a 

Z M (h) := f 7r M (d6)(l-26) h . 

It is immediate to observe that the expression inside the curly bracket is 
independent on i: 




where x(S) = J2aes% a - The a b° ve expression can be simplified if we define 
T as the set of the positions of the vector x different from zero. Indeed, the 
function x(S) can be expressed as 

x(S) = ]SnT[ 

where }Q[ = 1 if the cardinality of Q is odd and zero otherwise. Thus, 
replacing the sum over x with the sum over Y^Tc{k] m ■^ c l- El we fi nan y 
obtain 
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